// wave particle interaction

#include "toroidal.h"
#include <math.h>

double world::wpi_single(particle pa, double omega)
{
    double* eri_data = new double[Nx*Ny];
    double* ezi_data = new double[Nx*Ny];
    double* erq_data = new double[Nx*Ny];
    double* ezq_data = new double[Nx*Ny];
    this->partr_noiq(this->phies->get_fi(), eri_data);
    this->partz_noiq(this->phies->get_fi(), ezi_data);
    this->partr_noiq(this->phies->get_fq(), erq_data);
    this->partz_noiq(this->phies->get_fq(), ezq_data);
    field eri = field(eri_data);
    field ezi = field(ezi_data);
    field erq = field(erq_data);
    field ezq = field(ezq_data);

    int n = this->phies->n;
    double energy = 0.0;
    omega = omega/vh_vs_vA;
    double dt = 1.0/omega/100.0;
    int nrecord = 500;
    int intv = 10;
    double* record_R = new double[nrecord];
    double* record_Z = new double[nrecord];
    double* record_phi = new double[nrecord];
    double* record_time = new double[nrecord];
    double* pwr = new double[nrecord];
    // particle push code
    for(int i = 0; i < nrecord*intv; i++)
    {
        double vpar = pa.getvpara();
        cylv ccp = pa.getcrd();
        vec cp_xyz = cyl2xyz(ccp);
        pfvec em1 = this->emfield(ccp);
        phsp dX_phsp1 = this->push_rk2(em1, &pa, dt, vpar);
        cylv p1 = xyz2cyl(cp_xyz + dX_phsp1.X);
        pfvec em2 = this->emfield(p1);
        phsp dX_phsp2 = this->push_rk2(em2, &pa, dt, vpar+dX_phsp1.vp);
        cylv crdn = xyz2cyl(cp_xyz + 0.5*dX_phsp1.X + 0.5*dX_phsp2.X);
        pa.setcrd(crdn, vpar+0.5*dX_phsp1.vp+0.5*dX_phsp2.vp);
        if(!checkboundary(pa.getcrd()))
            printf("\ninside!!");
        if(i % intv == 0)
        {
            double tphase = 2.0+0.93*omega*dt*i;
            double cphase = 0.0 + 1.0*n*crdn.p;
            double er = eri.interp(crdn)*cos(tphase)*cos(cphase)+erq.interp(crdn)*sin(tphase)*sin(cphase);
            double ez = ezi.interp(crdn)*cos(tphase)*cos(cphase)+ezq.interp(crdn)*sin(tphase)*sin(cphase);
            double vr = rotate(0.5*(dX_phsp1.X + dX_phsp2.X), crdn.p).x;
            double power = er*vr + ez*0.5*(dX_phsp1.X + dX_phsp2.X).z;
            record_R[i/intv] = crdn.r;
            record_Z[i/intv] = crdn.z;
            pwr[i/intv] = power;
            record_phi[i/intv] = crdn.p;
            record_time[i/intv] = i*dt;
            energy += power;
        }
    }

    FILE* fp = fopen("/home/kun/Desktop/trace.csv","w");
    fprintf(fp,"R,Z,power,p, time\n");
    for(int i = 0; i < nrecord; i++)
    {
        fprintf(fp,"%f,%f,%f,%f,%f\n",record_R[i], 
        record_Z[i],pwr[i]*10000.0,record_phi[i],record_time[i]);
    }
    fclose(fp);
    printf("energy transfer:%f\n",1.0e4*energy); 
    return energy;
}

void world::wpi_multis_kc(double omega) // kinetic compression term
{
    double* eri_data = new double[Nx*Ny];
    double* epi_data = new double[Nx*Ny];
    double* ezi_data = new double[Nx*Ny];
    double* erq_data = new double[Nx*Ny];
    double* epq_data = new double[Nx*Ny];
    double* ezq_data = new double[Nx*Ny];
    this->partr_noiq(this->phies->get_fi(), eri_data);
    this->partz_noiq(this->phies->get_fi(), ezi_data);
    this->partr_noiq(this->phies->get_fq(), erq_data);
    this->partz_noiq(this->phies->get_fq(), ezq_data);
    vfield Ei = vfield(eri_data, epi_data, ezi_data);
    vfield Eq = vfield(erq_data, epq_data, ezq_data);
    vfiq Eperturb = vfiq(&Ei, &Eq, this->phies->n);
    vfiq* E0temp;

    double* Fr_data = new double[Nx*Ny];
    double* Fz_data = new double[Nx*Ny];
    partialx(*fast_distrib, Fr_data);
    partialy(*fast_distrib, Fz_data);
    field prf = field(Fr_data);
    field pzf = field(Fz_data);
    // prf.vtkoutput("/home/kun/Desktop/prF.vtk");
    
    int n = this->phies->n;
    double energy = 0.0;
    omega = omega/vh_vs_vA;
    double dt = 1.0/omega/100.0;
    int nrecord = 1;
    int intv = 10;
    double power = 0.0;
    // particle push code
    for(int i = 0; i < nrecord*intv; i++)
    {
        bool record = false;
        if(i % intv == 0)
        {
            power = 0.0;
            record = true;
        }
        for(int j = 0; j < this->fastptcs->npt; j++)
        {
            if(fastptcs->inside[j] == false)
                continue;
            particle pa = this->fastptcs->ppts[j];
            double vpar = pa.getvpara();
            cylv ccp = pa.getcrd();
            vec cp_xyz = cyl2xyz(ccp);
            pfvec em1 = this->emfield(ccp);
            phsp dX_phsp1 = this->push_rk2(em1, &pa, dt, vpar);
            cylv p1 = xyz2cyl(cp_xyz + dX_phsp1.X);
            pfvec em2 = this->emfield(p1);
            phsp dX_phsp2 = this->push_rk2(em2, &pa, dt, vpar+dX_phsp1.vp);
            cylv crdn = xyz2cyl(cp_xyz + 0.5*dX_phsp1.X + 0.5*dX_phsp2.X);
            pa.setcrd(crdn, vpar+0.5*dX_phsp1.vp+0.5*dX_phsp2.vp);
            fastptcs->inside[j] = checkboundary_core(crdn);
            
            E0temp = this->E0;
            this->E0 = &Eperturb;
            pfvec empt = this->emfield(ccp);
            phsp dX_pert = this->push_rk2(empt, &pa, dt, vpar);
            vec dz1 = rotate(dX_pert.X + -1.0*dX_phsp1.X, ccp.p);
            double floc = fast_distrib->interp(ccp);
            fastptcs->wdelta[j] += prf.interp(ccp)*dz1.x/floc;
            fastptcs->wdelta[j] += pzf.interp(ccp)*dz1.z/floc;
            this->E0 = E0temp;
        }
    }
} 

void world::kcp_single(particle pa, double omega)
{
    /* p0=(3.85,0.0,0.3) vb=0.4 mu=0.01
        omega= 0.03 t=115, tp=3.45
        n=3           */
    vfield E1in = field(phies->get_fi()).gradient();
    vfield E1qd = field(phies->get_fq()).gradient();
    this->E1i = &E1in;
    this->E1q = &E1qd;
    vfield gradF = this->fast_distrib->gradient();
    field Fr = field(gradF.getfr());
    field Fz = field(gradF.getfz());
    double wdelta = 0.0;
    double dt = 0.01;
    //int Nt = 1*11500;
    int Nt = 3000;
    int intv = 50;
    omega = 0.09;
    struct rcds
    {
        double rp;
        double zp;
        double tolp;
        double vEperp;
        double vEpara;
        double wdelta;
    };
    rcds* rcda = new rcds[Nt/intv];
    for(int i = 0; i < Nt; i++)
    {
        double tphase = omega*i*dt;
        double vpar = pa.getvpara();
        cylv ccp = pa.getcrd();
        vec cp_xyz = cyl2xyz(ccp);
        pfvec em0 = this->emfield(ccp);
        pfvec em1 = this->emfield_time(ccp, tphase);
        phsp dX_h0 = this->push_rk2(em0, &pa, dt, vpar);
        phsp dX_h1 = this->push_rk2(em1, &pa, dt, vpar);
        vec edrift = dX_h1.X + -1.0*dX_h0.X;
        cylv crdn = xyz2cyl(cp_xyz + dX_h0.X);
        pa.setcrd(crdn, vpar+dX_h0.vp);
        if(!checkboundary(pa.getcrd()))
        {    
            printf("\noutside!!");
            break;
        }

        double vEz = edrift.z;
        double vEr = edrift.x*cos(ccp.p) +edrift.y*sin(ccp.p);
        double fzp = fast_distrib->interp(ccp);
        double Frzp = Fr.interp(ccp);
        double Fzzp = Fz.interp(ccp);
        wdelta += (Frzp*vEr + Fzzp*vEz)/fzp;
        if (i % intv == 0)
        {
            double br = this->b->interp(ccp).x;
            double bz = this->b->interp(ccp).z;
            rcda[i/intv].tolp = 2*pi+fmod(tphase - phies->n*ccp.p,2.0*pi);
            rcda[i/intv].rp = ccp.r;
            rcda[i/intv].zp = ccp.z;
            
            double perpdrift = vEz*br - vEr*bz;
            double paradrift = vEr*br + vEz*bz;
            rcda[i/intv].vEpara = paradrift;
            rcda[i/intv].vEperp = perpdrift;
            rcda[i/intv].wdelta = wdelta;
        }
        // printf("deltaw:%f\n",vEz);
    }
    FILE* fp = fopen("/home/kun/Desktop/marker.csv","w");
    fprintf(fp,"R,Z,phase,vEperp,vEpara, dw\n");
    for(int i = 0; i < Nt/intv; i++)
    {
        fprintf(fp,"%f,%f,%f,%f,%f,%f\n",rcda[i].rp , rcda[i].zp,
        rcda[i].tolp, 3e5*rcda[i].vEperp, 3e5*rcda[i].vEpara,5*rcda[i].wdelta);
    }
    fclose(fp);
    cylv fnl = pa.getcrd();
    printf("\nfinal wdelta:%f,%f,%f,%f,%f\n",fnl.r,fnl.p,fnl.z,pa.getvpara(),wdelta);

}